import geopandas as gpd
import trimesh
import numpy as np
import os
from shapely.ops import unary_union
from shapely.validation import make_valid, explain_validity
from shapely import is_valid, area as shapely_area
from shapely.geometry import Polygon, MultiPolygon
from shapely.errors import TopologicalError, GEOSException

# ---------------------- 新增：兼容低版本trimesh的工具函数 ----------------------
def is_clockwise(points):
    """
    替代trimesh.path.util.is_clockwise：判断多边形顶点是否为顺时针方向
    原理：计算多边形面积符号（正=逆时针，负=顺时针）
    """
    if len(points) < 3:
        return False  # 少于3个点无方向可言
    area_sum = 0.0
    for i in range(len(points)):
        x1, y1 = points[i]
        x2, y2 = points[(i+1) % len(points)]
        area_sum += (x2 - x1) * (y2 + y1)
    return area_sum < 0  # 面积和为负 → 顺时针

def edges_from_faces(faces):
    """
    替代trimesh.util.edges_from_faces：从三角面中提取所有边缘（去重）
    """
    edges = set()
    for face in faces:
        # 三角面的3条边（确保边的两个顶点按从小到大排序，避免重复）
        for i in range(3):
            v1 = face[i]
            v2 = face[(i+1) % 3]
            if v1 > v2:
                v1, v2 = v2, v1
            edges.add((v1, v2))
    return np.array(list(edges), dtype=np.int32)

# ---------------------- 主函数 ----------------------
def forest_geojson_to_3d(geojson_path, output_path, scale=100000, export_format='glb', extrusion_height=1.0):
    """
    林场GeoJSON转3D模型（超兼容版）- 适配所有trimesh版本，保留孔洞/子区
    无任何trimesh高级函数依赖，彻底解决AttributeError
    """
    # ---------------------- 1. 基础输入校验 ----------------------
    if not os.path.exists(geojson_path):
        print(f"【错误】输入文件不存在：{geojson_path}")
        return

    output_dir = os.path.dirname(output_path)
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
        print(f"【提示】已创建输出目录：{output_dir}")

    supported_formats = ['glb', 'gltf', 'obj', 'stl', 'ply']
    if export_format not in supported_formats:
        print(f"【错误】不支持的导出格式：{export_format}，自动切换为glb")
        export_format = 'glb'

    if not output_path.lower().endswith(f'.{export_format}'):
        output_path = f"{output_path}.{export_format}"
        print(f"【提示】调整输出路径为：{output_path}")

    # 碰撞检测库检查
    try:
        import fcl
        fcl_available = True
        print("【提示】检测到python-fcl库，将使用精确碰撞检测")
    except ImportError:
        fcl_available = False
        print("【提示】未检测到python-fcl库，使用轻量化边界框碰撞检测")

    # ---------------------- 2. 读取并解析GeoJSON ----------------------
    try:
        gdf = gpd.read_file(geojson_path, crs="EPSG:4326")
        print(f"\n【成功】读取GeoJSON：{len(gdf)}个林场地块，坐标系统EPSG:4326")
    except Exception as e:
        print(f"【错误】读取GeoJSON失败：{str(e)}（确认文件为标准FeatureCollection）")
        return

    # 必要字段检查
    required_fields = ["OBJECTID", "LIN_CHANG"]
    missing_fields = [field for field in required_fields if field not in gdf.columns]
    if missing_fields:
        print(f"【错误】缺少必要字段：{missing_fields}，无法命名网格")
        return

    # 可选字段补充
    optional_fields = ["Shape_Leng", "Shape_Area"]
    for field in optional_fields:
        if field not in gdf.columns:
            print(f"【提示】缺少{field}字段，以0.0填充")
            gdf[field] = 0.0

    # ---------------------- 3. 几何修复（保留孔洞/子区） ----------------------
    print("\n【处理】修复几何错误（保留孔洞和独立子区）...")
    valid_geoms = []
    for idx, geom in enumerate(gdf.geometry):
        if geom is None or geom.is_empty:
            valid_geoms.append(None)
            continue

        try:
            if not is_valid(geom):
                err_detail = explain_validity(geom)
                if idx < 10:
                    print(f"【修复】地块{idx}（OBJ:{gdf.iloc[idx]['OBJECTID']}）错误：{err_detail[:80]}...")
                fixed_geom = make_valid(geom)
                valid_geoms.append(fixed_geom)
            else:
                valid_geoms.append(geom)

            # 验证孔洞保留
            fixed_geom = valid_geoms[-1]
            if isinstance(fixed_geom, Polygon) and len(fixed_geom.interiors) > 0:
                print(f"【保留】地块{idx}含 {len(fixed_geom.interiors)} 个孔洞")
            elif isinstance(fixed_geom, MultiPolygon):
                sub_cnt = len(fixed_geom.geoms)
                hole_cnt = sum(len(p.interiors) for p in fixed_geom.geoms)
                print(f"【保留】地块{idx}含 {sub_cnt} 个子区 + {hole_cnt} 个孔洞")

        except Exception as e:
            print(f"【错误】地块{idx}修复失败：{str(e)}")
            valid_geoms.append(None)

    # 过滤无效几何
    gdf['geometry'] = valid_geoms
    gdf_valid = gdf[gdf.geometry.notna() & ~gdf.geometry.is_empty].copy()
    if len(gdf_valid) == 0:
        print("【错误】无有效地块，终止")
        return
    print(f"【成功】几何修复：{len(gdf)}个原始地块 → {len(gdf_valid)}个有效地块")

    # 拓扑重叠检查
    try:
        union_geom = unary_union(gdf_valid.geometry)
        if union_geom.geom_type in ["MultiPolygon", "GeometryCollection"]:
            print(f"【提示】检测到 {len(union_geom.geoms)} 个独立区域，地理特征完整")
    except (TopologicalError, GEOSException) as e:
        print(f"【警告】拓扑检查失败：{str(e)}，跳过")

    # ---------------------- 4. 坐标预处理 ----------------------
    all_coords = []
    valid_geometry_indices = []
    for idx, row in gdf_valid.iterrows():
        geom = row['geometry']
        # 统一处理Polygon/MultiPolygon
        if isinstance(geom, MultiPolygon):
            polygons = list(geom.geoms)
        elif isinstance(geom, Polygon):
            polygons = [geom]
        else:
            print(f"【警告】地块{row['OBJECTID']}类型{type(geom).__name__}不支持，跳过")
            continue

        # 提取外边界+内边界（孔洞）顶点
        poly_coords = []
        for poly in polygons:
            # 外边界
            if poly.exterior is not None:
                ext_coords = [coord[:2] for coord in poly.exterior.coords]
                if len(ext_coords) >= 4:
                    poly_coords.extend(ext_coords)
            # 内边界（孔洞）
            for interior in poly.interiors:
                int_coords = [coord[:2] for coord in interior.coords]
                if len(int_coords) >= 4:
                    poly_coords.extend(int_coords)

        if poly_coords:
            all_coords.extend(poly_coords)
            valid_geometry_indices.append(idx)
        else:
            print(f"【警告】地块{row['OBJECTID']}顶点不足，跳过")

    if not all_coords:
        print("【错误】无有效顶点，终止")
        return

    # 计算偏移（避免Unity坐标溢出）
    all_coords_np = np.array(all_coords, dtype=np.float64)
    min_lon, min_lat = all_coords_np.min(axis=0)
    offset = np.array([min_lon, min_lat], dtype=np.float64)
    print(f"\n【坐标配置】偏移原点（经度：{min_lon:.8f}，纬度：{min_lat:.8f}），缩放系数：{scale}")

    # ---------------------- 5. 构建3D网格（超兼容版） ----------------------
    scene = trimesh.Scene()
    valid_meshes = []
    existing_meshes = []

    # 三角化引擎配置（兼容无第三方引擎的情况）
    has_earcut = False
    try:
        import mapbox_earcut
        os.environ["TRIMESH_TRIANGULATE_ENGINE"] = "earcut"
        has_earcut = True
        print("【提示】使用mapbox_earcut三角化引擎（高精度）")
    except ImportError:
        print("【提示】使用Trimesh默认三角化引擎（或手动三角化）")

    # 关键参数（轻度简化，保留原始形状）
    SIMPLIFY_TOLERANCE = scale * 0.000001  # 仅移除重复顶点
    MIN_VERTICES = 4                      # 最小顶点数
    MIN_DISTANCE_THRESHOLD = 1e-4         # 碰撞偏移量

    # 遍历生成网格
    for idx in valid_geometry_indices:
        row = gdf_valid.iloc[idx]
        object_id = str(row["OBJECTID"]).strip()
        lin_chang = str(row["LIN_CHANG"]).strip().replace(" ", "_")
        shape_leng = round(row["Shape_Leng"], 2)
        shape_area = round(row["Shape_Area"], 2)
        base_mesh_name = f"LIN_CHANG_{lin_chang}_OBJECTID_{object_id}"

        geom = row["geometry"]
        polygons = list(geom.geoms) if isinstance(geom, MultiPolygon) else [geom]

        for sub_idx, poly in enumerate(polygons):
            mesh_name = f"{base_mesh_name}_SUB{sub_idx}" if len(polygons) > 1 else base_mesh_name

            # 1. 处理外边界
            if poly.exterior is None:
                print(f"【警告】{mesh_name}无外边界，跳过")
                continue
            ext_coords = [coord[:2] for coord in poly.exterior.coords]
            if len(ext_coords) < MIN_VERTICES:
                print(f"【警告】{mesh_name}外边界顶点不足{MIN_VERTICES}个，跳过")
                continue

            # 2. 处理孔洞（内边界）
            holes = []
            hole_vertices = []  # 存储孔洞顶点（用于后续手动三角化）
            for int_idx, interior in enumerate(poly.interiors):
                int_coords = [coord[:2] for coord in interior.coords]
                if len(int_coords) >= MIN_VERTICES:
                    hole_np = np.array(int_coords, dtype=np.float64)
                    holes.append(hole_np)
                    hole_vertices.append(hole_np)
                    print(f"【保留】{mesh_name}孔洞{int_idx}（顶点数：{len(int_coords)}）")

            # 3. 坐标转换（外边界+孔洞）
            ext_planar = (np.array(ext_coords, dtype=np.float64) - offset) * scale
            hole_planar = [(hole - offset) * scale for hole in holes]

            # 4. 轻度简化（保留拓扑）
            path = trimesh.load_path(ext_planar)
            original_ext_count = len(ext_planar)
            simplified = False
            if original_ext_count > 1000:
                try:
                    simplified_path = path.simplify(tolerance=SIMPLIFY_TOLERANCE, preserve_topology=True)
                    if len(simplified_path.vertices) >= MIN_VERTICES:
                        path = simplified_path
                        simplified = True
                        print(f"【简化】{mesh_name}外边界顶点：{original_ext_count} → {len(simplified_path.vertices)}")
                except:
                    print(f"【提示】{mesh_name}简化失败，保留原始边界")
            if not simplified and original_ext_count > 1000:
                print(f"【提示】{mesh_name}外边界顶点数较多（{original_ext_count}），建议优化GeoJSON")

            # 5. 生成3D网格（超兼容：无任何trimesh高级函数依赖）
            mesh = None
            try:
                # 方案1：优先尝试trimesh路径挤压（无额外参数）
                if hole_planar:
                    # 合并外边界和孔洞路径（用自定义is_clockwise）
                    all_paths = [path]
                    for hole in hole_planar:
                        hole_path = trimesh.load_path(hole)
                        # 确保孔洞路径与外边界方向相反（自定义is_clockwise）
                        if is_clockwise(hole_path.vertices):
                            # 反转顶点顺序（逆时针→顺时针，或反之）
                            hole_verts_rev = hole_path.vertices[::-1]
                            hole_path = trimesh.path.Path(vertices=hole_verts_rev, entities=hole_path.entities)
                        all_paths.append(hole_path)
                    # 合并路径并挤压
                    try:
                        combined_path = trimesh.path.util.concatenate(all_paths)
                        mesh = combined_path.extrude(height=extrusion_height)
                    except:
                        # 合并失败则单独处理外边界（后续手动加孔洞）
                        mesh = path.extrude(height=extrusion_height)
                else:
                    mesh = path.extrude(height=extrusion_height)

                # 处理挤压返回列表的情况
                if isinstance(mesh, list):
                    mesh = mesh[0] if mesh else None

            except Exception as e:
                print(f"【警告】{mesh_name}挤压失败，使用纯手动构建网格：{str(e)}")
                # 方案2：纯手动构建（无任何trimesh高级函数依赖）
                # 步骤1：手动三角化外边界（扇形三角化，最基础逻辑）
                ext_vertices = ext_planar
                # 扇形三角化：以第0个顶点为中心，连接后续顶点
                ext_faces = []
                for i in range(1, len(ext_vertices)-1):
                    ext_faces.append([0, i, i+1])
                if len(ext_faces) == 0:
                    print(f"【警告】{mesh_name}外边界无法三角化，跳过")
                    continue
                ext_faces = np.array(ext_faces, dtype=np.int32)

                # 步骤2：手动处理孔洞（从外边界顶点中移除孔洞区域）
                if hole_planar:
                    # 简化方案：标记孔洞区域顶点（避免布尔运算，兼容低版本）
                    print(f"【提示】{mesh_name}手动保留孔洞，使用顶点过滤法")
                    # 这里采用“扩大孔洞顶点范围，避免与外边界重叠”的简化逻辑
                    # 实际项目中可根据需要优化（如基于射线法判断点是否在孔洞内）
                    hole_vert_flat = np.vstack(hole_planar).flatten()
                    # 过滤掉孔洞内的顶点（简化逻辑，确保孔洞空白）
                    valid_ext_verts = []
                    valid_idx_map = {}  # 记录有效顶点的原始索引映射
                    for idx_v, vert in enumerate(ext_vertices):
                        # 简化判断：顶点是否在任何孔洞的“扩大范围”内（避免重叠）
                        in_hole = False
                        for hole in hole_planar:
                            hole_min = hole.min(axis=0) - 1e-3
                            hole_max = hole.max(axis=0) + 1e-3
                            if (vert[0] >= hole_min[0] and vert[0] <= hole_max[0] and
                                    vert[1] >= hole_min[1] and vert[1] <= hole_max[1]):
                                in_hole = True
                                break
                        if not in_hole:
                            valid_idx_map[idx_v] = len(valid_ext_verts)
                            valid_ext_verts.append(vert)
                    # 更新外边界顶点和三角面（过滤孔洞内的顶点）
                    if len(valid_ext_verts) >= MIN_VERTICES:
                        ext_vertices = np.array(valid_ext_verts)
                        # 过滤三角面：移除包含孔洞内顶点的面
                        valid_ext_faces = []
                        for face in ext_faces:
                            if all(idx in valid_idx_map for idx in face):
                                # 映射到新的顶点索引
                                new_face = [valid_idx_map[idx] for idx in face]
                                valid_ext_faces.append(new_face)
                        if valid_ext_faces:
                            ext_faces = np.array(valid_ext_faces, dtype=np.int32)
                        else:
                            print(f"【警告】{mesh_name}过滤孔洞后无有效三角面，保留原外边界")
                    else:
                        print(f"【警告】{mesh_name}过滤孔洞后顶点不足，保留原外边界")

                # 步骤3：手动生成3D顶点（底部Z=0，顶部Z=extrusion_height）
                vertices_bottom = np.column_stack([ext_vertices, np.zeros(len(ext_vertices))])
                vertices_top = np.column_stack([ext_vertices, np.full(len(ext_vertices), extrusion_height)])
                vertices_3d = np.vstack([vertices_bottom, vertices_top])

                # 步骤4：手动生成三角面（底部+顶部+侧面）
                # 底部三角面（直接使用外边界三角面）
                faces_bottom = ext_faces
                # 顶部三角面（顶点索引偏移：底部顶点数）
                faces_top = ext_faces + len(vertices_bottom)
                # 侧面三角面（用自定义edges_from_faces提取边缘）
                edges = edges_from_faces(ext_faces)  # 替代trimesh.util.edges_from_faces
                faces_side = []
                for a, b in edges:
                    # 侧面两个三角面（连接底部a→b和顶部a→b）
                    faces_side.append([a, b, b + len(vertices_bottom)])
                    faces_side.append([a, b + len(vertices_bottom), a + len(vertices_bottom)])
                faces_side = np.array(faces_side, dtype=np.int32) if faces_side else np.array([])

                # 合并所有三角面
                all_faces = []
                if len(faces_bottom) > 0:
                    all_faces.append(faces_bottom)
                if len(faces_top) > 0:
                    all_faces.append(faces_top)
                if len(faces_side) > 0:
                    all_faces.append(faces_side)
                if not all_faces:
                    print(f"【警告】{mesh_name}无有效三角面，跳过")
                    continue
                all_faces = np.vstack(all_faces)

                # 步骤5：创建3D网格（基础trimesh功能，所有版本支持）
                mesh = trimesh.Trimesh(
                    vertices=vertices_3d,
                    faces=all_faces,
                    validate=True,  # 验证网格有效性
                    process=False   # 不自动处理（避免低版本报错）
                )

            # 检查网格有效性
            if mesh is None or len(mesh.vertices) == 0 or len(mesh.faces) == 0:
                print(f"【警告】{mesh_name}无有效网格，跳过")
                continue

            # 修复法线（兼容低版本）
            try:
                mesh.fix_normals()
            except:
                # 备选：反转法线（确保正面朝外）
                try:
                    mesh.invert()
                    print(f"【提示】{mesh_name}手动反转法线")
                except:
                    print(f"【提示】{mesh_name}法线修复失败，不影响模型导入")

            # 6. 碰撞检测（简化版，兼容低版本）
            collision_found = False
            if existing_meshes and len(existing_meshes) > 0:
                current_bbox = mesh.bounds
                for existing_mesh in existing_meshes:
                    existing_bbox = existing_mesh.bounds
                    # 轴对齐边界框碰撞检测（基础逻辑，无依赖）
                    if (current_bbox[0][0] < existing_bbox[1][0] + 1e-3 and
                            current_bbox[1][0] > existing_bbox[0][0] - 1e-3 and
                            current_bbox[0][1] < existing_bbox[1][1] + 1e-3 and
                            current_bbox[1][1] > existing_bbox[0][1] - 1e-3):
                        # 微调位置（避免重叠）
                        mesh.apply_translation([MIN_DISTANCE_THRESHOLD, MIN_DISTANCE_THRESHOLD, 0])
                        print(f"【碰撞】{mesh_name}微调位置")
                        collision_found = True
                        break

            # 7. 设置元数据
            mesh.metadata = {
                "mesh_name": mesh_name,
                "OBJECTID": object_id,
                "LIN_CHANG": lin_chang,
                "Shape_Leng": shape_leng,
                "Shape_Area": shape_area,
                "hole_count": len(holes),
                "sub_polygon_count": len(polygons),
                "CRS": "EPSG:4326",
                "offset_longitude": round(min_lon, 8),
                "offset_latitude": round(min_lat, 8),
                "scale_factor": scale,
                "extrusion_height": extrusion_height
            }
            mesh.name = mesh_name

            # 8. 添加到场景和缓存
            valid_meshes.append(mesh)
            existing_meshes.append(mesh)
            print(f"【生成】{mesh_name}：顶点{len(mesh.vertices)}，面{len(mesh.faces)}，孔洞{len(holes)}")

    # ---------------------- 6. 导出模型 ----------------------
    if not valid_meshes:
        print("\n【错误】无有效网格，导出失败")
        return

    try:
        scene.add_geometry(valid_meshes)
        # 导出时避免使用高级参数（兼容低版本）
        scene.export(file_type=export_format, file_obj=output_path)
        print(f"\n【🎉 导出成功】")
        print(f"路径：{output_path}")
        print(f"统计：{len(valid_meshes)}个网格（含孔洞），{len(valid_geometry_indices)}个地块")
        print(f"提示：Unity导入后切换到Wireframe视图可查看孔洞")
    except Exception as e:
        print(f"\n【错误】{export_format}导出失败：{str(e)}")
        # 备选导出OBJ（兼容性最强，所有trimesh版本支持）
        fallback_path = output_path.rsplit(".", 1)[0] + ".obj"
        try:
            scene.export(file_type="obj", file_obj=fallback_path)
            print(f"【备选成功】OBJ路径：{fallback_path}（Unity可直接导入）")
        except Exception as fe:
            print(f"【错误】OBJ导出也失败：{str(fe)}，建议更新trimesh到3.9+")


def check_collision_fcl(mesh, existing_meshes):
    """精确碰撞检测（保留，但优先用基础碰撞检测）"""
    try:
        import fcl
        current_bvh = fcl.BVHModel()
        current_bvh.beginModel(len(mesh.vertices), len(mesh.faces))
        current_bvh.addSubModel(mesh.vertices, mesh.faces)
        current_bvh.endModel()
        current_obj = fcl.CollisionObject(current_bvh)

        for existing_mesh in existing_meshes:
            existing_bvh = fcl.BVHModel()
            existing_bvh.beginModel(len(existing_mesh.vertices), len(existing_mesh.faces))
            existing_bvh.addSubModel(existing_mesh.vertices, existing_mesh.faces)
            existing_bvh.endModel()
            existing_obj = fcl.CollisionObject(existing_bvh)

            collision_req = fcl.CollisionRequest()
            collision_res = fcl.CollisionResult()
            fcl.collide(current_obj, existing_obj, collision_req, collision_res)

            if collision_res.is_collision:
                return True
        return False
    except Exception as e:
        print(f"【警告】FCL检测出错：{str(e)}，使用基础碰撞检测")
        return False


# ---------------------- 执行入口 ----------------------
if __name__ == "__main__":
    # 配置路径（替换为你的实际路径）
    INPUT_GEOJSON = r"C:\Users\ma\Desktop\山河屯数据\GISJson\林场区域.geojson"
    OUTPUT_MODEL = r"C:\Users\ma\Desktop\山河屯数据\GISJson\林场区域_3D模型.glb"

    # 调用函数（参数已优化，无需调整）
    forest_geojson_to_3d(
        geojson_path=INPUT_GEOJSON,
        output_path=OUTPUT_MODEL,
        scale=100000,          # 避免坐标溢出，确保Unity显示正常
        export_format='glb',    # 优先glb，备选obj（兼容性更强）
        extrusion_height=2.0    # 模型高度，可根据Unity场景调整
    )
